15  Gauss Quadrature

Note. These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Chapter 1: How computers add

Chapter 2: Solving nonlinear equations in 1d

Chapter 3: Polynomial interpolation

Chapter 4: Numerical Integration

Next - Chapter 5: Linear systems of equations

include("preamble.jl");
✓ file included! 

using: Plots, LaTeXStrings, Polynomials, PrettyTables 

Functions included: 
    simple_iteration, 
    Newton, 
    orderOfConvergence, 
    ChebyshevNodes 

Use @doc <<function>> for help

15.1 Previously…

Recall that we are considering quadrature rules of the following form:

\begin{align} \int_{a}^b f &\approx \sum_j w_j f(x_j). \end{align}

where \{w_j\} are the quadrature weights, \{x_j\} are the quadrature nodes, and (\{w_j\}, \{x_j\}) or \sum_{j} w_j f(x_j) is a quadrature rule.

Suppose p is the polynomial of degree \leq n interpolating f on X = \{x_0,\dots,x_n \}: that is p(x) = \sum_{j=0}^n \ell_j(x) f(x_j). We simply approximate the integral of f by instead integrating p:

\begin{align} \int_{a}^b f &\approx \int_a^b p = \sum_{j=0}^n \left[\int_a^b\ell_j(x)\mathrm{d}x \right] f(x_j) % \tag{$\star$} \end{align}

which gives a quadrature rule with weights w_j := \int_a^b\ell_j(x)\mathrm{d}x and nodes x_j.

15.1.1 Newton-Cotes

When the interpolation nodes \{x_0,\dots,x_n\} are equispaced, this is a Newton-Cotes quadrature rule.

We saw the following examples:

  • Rectangular rule:

\begin{align} \int_a^b f \approx (b-a) f(a) \end{align}

  • Midpoint rule (from exam):

\begin{align} \int_a^b f \approx (b-a) f\big( \tfrac{a+b}2 \big) \end{align}

  • Trapezoid rule:

\begin{align} \int_a^b f \approx \frac{b-a}{2} \left( f(a) + f(b) \right) \end{align}

  • Simpson’s rule:

\begin{align} \int_a^b f \approx \frac{b-a}{6} \left( f(a) + 4 f\big( \tfrac{a+b}2 \big) + f(b) \right) \end{align}

Definition.

We say the quadrature rule (\{w_j\}, \{x_j\}) is exact for all polynomials of degree N if

\begin{align} \int_a^b p(x) \mathrm{d}x = \sum_{j=0}^n w_j p(x_j) \end{align}

for all p \in \mathcal P_N (all polynomials of degree \leq N).

Let us consider the quadrature rule (\star). For all polynomials P of degree \leq n, we have

\begin{align} \sum_{j=0}^n w_j P(x_j) &= \sum_{j=0}^n \left[ \int_a^b \ell_j(x) \mathrm{d}x\right] P(x_j) \\ % &= \int_a^b \sum_{j=0}^n \ell_j(x) P(x_j) \mathrm{d}x \\ % &= \int_a^b P \end{align}

In the final line, we have used the fact that P(x) = \sum_{j=0}^n \ell_j(x) P(x_j) beacuse P has degree \leq n. That is, the Newton-Cotes quadrature rules given by (\star) are exact for all polynomials of degree n.

Therefore:

  • Rectangular rule (n=0) is exact for all polynomials of degree 0
  • Trapezoid rule (n=1) is exact for all polynomials of degree 1

You have also seen that:

  • Midpoint rule (n=0) is exact for all polynomials of degree 1
  • Simpson rule (n=2) is exact for all polynomials of degree 3

This is a general fact: Suppose n is even and X = \{x_0, \dots, x_n\} is a set of interpolation nodes symmetric about \frac{a+b}{2} (that is, x_j + x_{n-j} = a+b ). Then, the quadrature rule (\star) is exact for all polynomials of degree n+1.

Suppose we have the quadrature rule (\star) with equispaced points X = \{x_0,\dots,x_n\} and it is exact for all polynomials of degree N (N = n for odd n and N = n+1 for even n). Let P the polynomial interpolant of f on \tilde{X} such that P has degree N and \tilde{X} contains X (i.e. P(x_j) = f(x_j) for all j=0,\dots,n). Then,

\begin{align*} \left| \int_a^b f - \sum_{j=0}^n w_j f(x_j) \right| % &= \left| \int_a^b (f - P ) \right| \\ % &\leq \frac{\|f^{(N+1)}\|_{L^\infty([a,b])}}{(N+1)!} \int_a^b \left| \ell_{\tilde{X}}(x) \right| \mathrm{d}x \end{align*}

Composite Rules. Fix K and define h := \frac{b-a}{K}. Then, split [a,b] into the K intervals [a + kh, a + (k+1)h] for k = 0,\dots,K-1:

\begin{align} \int_a^b f = \sum_{ k=0 }^{K-1} \int_{a + kh}^{a + (k+1)h} f \end{align}

and apply the quadrature rule on each subinterval [a + kh, a + (k+1)h]. Using the error estimate on each sub-interval and summing over k, we saw that errors in the composite rules are as follows:

  • Rectangular: O(h),
  • Midpoint, Trapezoid: O(h^2),
  • Simpson: O(h^4)

15.2 L18-19: Gaussian Quadrature

Remark.

By considering the change of variables y = \frac{2}{b-a} \left( x - \frac{a + b}{2} \right) and defining \tilde{f}(y) = \frac{b-a}{2} f\big( \tfrac{b-a}2 y + \tfrac{a+b}2 \big), we have

\begin{align} \int_a^b f(x) \mathrm{d}x = \int_{-1}^1 \tilde{f}(y) \mathrm{d}y. \end{align}

Moreover, if p is the polynomial interpolation of f on X \subset [a,b] then \tilde{p} is the polynomial interpolation of \tilde{f} on \tilde{X} \subset [-1,1] where \tilde{X} = \{ \frac{2}{b-a} \left( x - \frac{a + b}{2} \right) : x \in X \}. As a result, we may without loss of generality consider the interval [-1,1].

Before: Fix X find \{w_j\} such that the quadrature rule

\begin{align} \int_a^b f \approx \sum_{j} w_j f(x_j) \end{align}

is exact on \mathcal P_n for maximal n.

Idea: Choose X and \{w_j\}.

Example.

Choose w_0, w_1 and x_0, x_1 such that

\begin{align} \int_{-1}^1 f \approx w_0 f(x_0) + w_1 f(x_1) \end{align}

is exact for all polynomials of degree 3 = 2n+1.

Answer.

Using the functions, f = 1, x, x^2, x^3, we want to solve

\begin{align} 2 &= w_0 + w_1 \\ 0 &= w_0 x_0 + w_1 x_1 \\ \frac{2}{3} &= w_0 x_0^2 + w_1 x_1^2 \\ 0 &= w_0 x_0^3 + w_1 x_1^3 \end{align}

We have

\begin{align} 0 &= ( w_0 x_0 + w_1 x_1 )( x_0 + x_1 ) \\ &= (w_0 x_0^2 + w_1 x_1^2) + ( w_0 + w_1 ) x_0 x_1 \\ &= \frac23 + 2 x_0 x_1 \end{align}

and so x_0 x_1 = -\frac13. Similarly, we have

\begin{align} \frac23 (x_0 + x_1) &= ( w_0 x_0^2 + w_1 x_1^2 )( x_0 + x_1 ) \\ &= (w_0 x_0^3 + w_1 x_1^3) + ( w_0 x_0 + w_1 x_1 ) x_0 x_1 \\ &= 0 + 0 x_0 x_1 = 0 \end{align}

and so x_0 = - x_1 and x_0 = - \frac{1}{\sqrt{3}}, x_1 = \frac{1}{\sqrt{3}}.

Finally, we have 0 = w_0 x_0 + w_1 x_1 = (w_0 - w_1) x_0 and so w_0 = w_1 = 1.

The above example gives the quadrature rule

\begin{align} \int_{-1}^1 f \approx f\big( -\tfrac{\sqrt3}{3} \big) + f\big( \tfrac{\sqrt3}{3} \big) \end{align}

which is exact for all polynomials of degree \leq 3.

Q: How to generalise this??

It turns out that one can choose the interpolation nodes and weights so that the quadrature rule

\begin{align} \int_{-1}^1 f \approx \sum_{j=0}^n w_j f(x_j) \end{align}

is exact for all polynomials of degree 2n+1. The solution in the general case requires a set of orthogonal polynomials:

15.3 Legendre Polynomials

Definition.

Let the Legendre polynomial P_n be the monic polynomial of degree n such that

\begin{align} (P_n, q)_{L^2} := \int_{-1}^1 P_n(x) q(x) \mathrm{d}x = 0 \end{align}

for all q \in P_{n-1}.

Remark.

We say that f is orthogonal to g if (f,g)_{L^2} = 0.

Remark (An equivalent definition).

P_n is monic and (P_n, P_m) = 0 for all n \not= m.

It may not be immediately clear that such polynomials exist, but we may compute the first few explicitly:

Example.

When n=0, there is only one monic polynomial: P_0(x) = 1,

When n=1, we require P_1(x) = x + a to satisfy

\begin{align} 0 = \int_{-1}^1 P_1(x) 1 \mathrm{d}x = 2 a \end{align}

and so P_1(x) = x.

When n=2, we have P_2(x) = x^2 + a x + b and so

\begin{align} 0 &= \int_{-1}^1 P_2(x) 1 \mathrm{d}x = 2b + \frac23 \\ 0 &= \int_{-1}^1 P_2(x) x \mathrm{d}x = \frac{2a}3 \end{align}

therefore P_2(x) = x^2 - \frac13.

Theorem.

The Legendre polynomials exist for all n.

(we will come back to the proof)

Remark.

The roots of P_2 are \{ -\tfrac{\sqrt3}{3}, \tfrac{\sqrt3}{3} \} (which are the nodes used in the quadrature rule we had above)

Recall P_n are the Legendre polynomials and define the set of Legendre points X by

\begin{align} X := \{ \text{zeros of } P_{n+1} \}. \end{align}

Claim.

  • X is a set of n+1 distinct points in [-1,1]

Proof.

Notice that P_{n+1}(x) = \prod_{j=0}^n(x-x_j) (i.e. the monic polynomial with roots X).

Suppose that x_0 = x_1 and define q(x) := \prod_{j=2}^n (x - x_j). Then, since P_{n+1} is orthogonal to q (as q is of degree n-1), we have

\begin{align} 0 &= \int_{-1}^1 P_{n+1}(x) q(x) \mathrm{d}x \\ % &= \int_{-1}^1 \left(\prod_{j=0}^n (x - x_j) \right) \left( \prod_{j=2}^n (x - x_j) \right) \mathrm{d}x \\ % &= \int_{-1}^1 (x - x_0)^2 \prod_{j=2}^n (x - x_j)^2 \mathrm{d}x >0, \end{align}

a contradiction (the integral of a non-negative polynomial on a positive interval is positive). By relabeling the indices in X = \{x_0,\dots,x_n\}, we see that X is a set of pairwise distinct nodes.

Suppose x_k \not\in [-1,1]. Then, on defining q(x) := \prod_{j=0 \,:\, j \not= k }^n (x - x_j), a polynomial of degree n, we have

\begin{align} 0 &= \int_{-1}^1 P_{n+1}(x) q(x) \mathrm{d}x \\ % &= \int_{-1}^1 \left(\prod_{j=0}^n (x - x_j) \right) \left( \prod_{j\not=k} (x - x_j) \right) \mathrm{d}x \\ % &= \int_{-1}^1 (x - x_k) \prod_{j\not=k} (x - x_j)^2 \mathrm{d}x. \end{align}

Now, since x_k \not\in [-1,1], we have (x - x_k) is either strictly positive or negative for x \in [-1,1]. Therefore, \int_{-1}^1 (x - x_k) \prod_{j\not=k} (x - x_j)^2\mathrm{d}x is either strictly positive or negative, a contradiction. As a result, each x_k is in the interval [-1,1].

Definition (Gaussian Quadrature)

Let p be the polynomial interpolation of f on the n+1 Legendre points. Then, we may define the quadrature rule

\begin{align} \int_{-1}^1 f \approx \int_{-1}^{1} p = \sum_{j=0}^n w_j f(x_j) \end{align}

Theorem.

Gaussian quadrature is exact for polynomials of degree 2n+1.

Proof.

First, note that we have already seen that Gaussian quadrature is exact for all polynomials of degree n.

Let X = \{x_0,\dots,x_n\} be the zeros of P_{n+1}, w_j the corresponding quadrature nodes, and let P be a polynomial of degree 2n+1. On dividing P by P_{n+1}, there exists q_n, r_n \in \mathcal P_n such that

\begin{align} P(x) = P_{n+1}(x) q_n(x) + r_n(x). \end{align}

Notice that P(x_j) = P_{n+1}(x_j) q_n(x_j) + r_n(x_j) = r_n(x_j) so that r_n is a degree n polynomial interpolant of P on X. Using this, combined with the fact that P_{n+1} is orthogonal to q_n, we have

\begin{align} \int_{-1}^1 P &= \int_{-1}^{1} \left[ P_{n+1}(x) q_n(x) + r_n(x) \right] \mathrm{d}x \\ &= 0 + \int_{-1}^{1} r_n(x) \mathrm{d}x \\ % &= \sum_{j=0}^n w_j r_n(x_j) \\ % &= \sum_{j=0}^n w_j P(x_j) % \end{align}

Therefore, Gaussian quadrature is exact for all polynomials of degree \leq 2n+1.

Claim.

  • \sum_{j=0}^n w_j = 2
  • w_j \geq 0,

Proof.

Notice that \sum_{k=0}^n \ell_k(x) is a polynomial of degree n which is equal to 1 on the set of n+1 points X and so \sum_{k=0}^n \ell_k(x)= 1 for all x. Therefore, because the quadrature rule is exact for polynomials of degree 2n+1, we have

\begin{align} 2 &= \int_{-1}^1 \mathrm{d}x \\ &= \int_{-1}^1 \sum_{k=0}^n \ell_k(x) \mathrm{d}x \\ &= \sum_{j=0}^n w_j \left[ \sum_{k=0}^n \ell_j(x_k) \right] \\ &= \sum_{j=0}^n w_j. \end{align}

Moreover, \ell_k(x)^2 is a polynomial of degree 2n and the quadrature rule is exact for all polynomials of degree 2n +1, we have

\begin{align} 0 &\leq \int_{-1}^1 \ell_k(x)^2 \mathrm{d}x \\ &= \sum_j w_j \left[ \ell_k(x_j)^2 \right] \\ &= w_j \end{align}

Therefore, the weights are positive and sum to 2 = \int_{-1}^1 1 \mathrm{d}x.

Theorem.

The Legendre polynomials P_n exist.

Proof. (Presentation)

We will use the notation

\begin{align} \|f\|_{L^2} := \left( \int_{-1}^1 |f(x)|^2 \mathrm{d}x \right)^{1/2}. \end{align}

Let p_0 = 1 and p_1 = x.

If p_0, p_1, \dots, p_n have been defined, we define

\begin{align} a_n &:= \frac{ \int_{-1}^1 x p_{n}(x)^2 \mathrm{d}x }{\|p_{n}\|_{L^2}^2} \qquad \text{and}\\ b_{n} &:= \frac{ \int_{-1}^1 x p_{n}(x) p_{n-1}(x) \mathrm{d}x }{\|p_{n-1}\|_{L^2}^2} \end{align}

and thus p_{n+1}(x) := (x - a_n) p_n(x) - b_n p_{n-1}(x) for all n \geq 1.

We claim that p_n(x) is the n^\text{th} Legendre polynomial. We prove this by induction: the statement is true for n=0 and n=1. Suppose the statement is true for n.

Then, notice that p_{n+1}(x) is a monic polynomial of degree n+1.

Moreover, we have: for all j < n-1,

\begin{align} \int_{-1}^1 p_{n+1}(x) p_{j}(x) \mathrm{d}x % &= \int_{-1}^1 \left[ (x - a_n) p_{n}(x) - b_n p_{n-1} \right] p_j(x) \mathrm{d}x \\ % &= \int_{-1}^1 p_{n}(x) \left[ (x - a_n) p_j(x) \right] - b_n \int_{-1}^1 p_{n-1} p_j \\ % &= 0. \end{align}

Here, we have used the fact that (x-a_n) p_j is a polynomial of degree j+1 < n, p_j is a polynomial of degree j < n-1, and p_n is the n^\text{th} Legendre polynomial.

Moreover, we have

\begin{align} \int_{-1}^1 p_{n+1}(x) p_{n-1}(x) \mathrm{d}x % &= \int_{-1}^1 \left[ (x - a_n) p_{n}(x) - b_n p_{n-1} \right] p_{n-1}(x) \mathrm{d}x \\ % &= \int_{-1}^1 x p_{n}(x) p_{n-1}(x) - b_n \int_{-1}^1 | p_{n-1}(x) |^2 \mathrm{d}x \\ % &= 0 \end{align}

and

\begin{align} \int_{-1}^1 p_{n+1}(x) p_{n}(x) \mathrm{d}x % &= \int_{-1}^1 \left[ (x - a_n) p_{n}(x) - b_n p_{n-1} \right] p_{n}(x) \mathrm{d}x \\ % &= \int_{-1}^1 (x - a_n) p_{n}(x) p_{n}(x) \\ % &= 0 \end{align}

As a result, we have p_{n+1} is monic and orthogonal to all polynomials of degree \leq n.

We may use the recursion in the proof to compute the Legendre polynomials: here we plot P_n for n \leq 8 and the set of Legendre points

Next we approximate the integral

\begin{align} \int_{-1}^{+1} \frac{1}{1 + e^{10(x - 0.5)}} \mathrm{d}x \end{align}

by using (i) the composite trapezoid rule, (ii) composite simpson rule (from last time) and (iii) using Newton quadrature.

15.4 Error estimate for n=1

We have X = \{x_0, x_1\} = \{ -\frac{\sqrt3}{3}, \frac{\sqrt3}3 \}. Quadrature in these points is exact for all polynomials of degree 2(1) + 1 = 3. We therefore let p be the polynomial of degree 3 such that p(x_j) = f(x_j) and p'(x_j) = f'(x_j) for j=0,1 and notice that the error is

\begin{align} &\left| \int_{-1}^1 f - \left( f(x_0) + f(x_1) \right) \right| \\ % &= \left| \int_{-1}^1 f - \left( p(x_0) + p(x_1) \right) \right| \\ % &= \left| \int_{-1}^1 (f - p) \right|. \end{align}

Since p is the Hermite polynomial interpolant on \{\{ x_0, x_0, x_1, x_1 \}\}, there exists \xi_x such that f(x) - p(x) = \frac{f^{(4)}(\xi_x)}{4!} (x-x_0)^2 (x-x_1)^2 [this is the generalisation of the Taylor remainder theorem we saw when we were considering Hermite interpolation]. Therefore, using the mean value theorem, there exists \xi \in [-1,1] such that

\begin{align} &\left| \int_{-1}^1 (f - p) \right| \\ &= \left| \frac{f^{(4)}(\xi)}{4!} \right| \int_{-1}^1 (x-x_0)^2 (x-x_1)^2 \mathrm{d}x \\ &= \frac{8}{45} \frac{1}{4!} \left| f^{(4)}(\xi) \right| \\ &= \frac{1}{135} \left| f^{(4)}(\xi) \right|. \end{align}

Example.

Approximate \log 2 = \int_1^2 \frac{1}{x} \mathrm{d}x using Gauss quadrature and evaluate an upper bound for the error.

Answer.

First notice that

\begin{align} \log 2 &= \int_1^2 \frac1x \mathrm{d}x \\ % &= \frac{1}{2} \int_{-1}^1 \frac1{ \frac{1}2 x + \frac32 } \mathrm{d}x \\ % &= \int_{-1}^1 \frac1{ 3 + x } \mathrm{d}x \\ % &\approx \frac1{ 3 - \frac{\sqrt{3}}3 } + \frac1{ 3 + \frac{\sqrt{3}}3 } \\ % &= \frac9{13} \end{align}

The error here is \log 2 - \tfrac{9}{13} \approx 8.39 \times 10^{-4}.

The error is \frac{1}{135} f^{(4)}(\xi) where f(x) = \frac{1}{3+x} and \xi is some point in [-1,1]. Taking derivatives of f, we find that f^{(4)}(x) = \frac{24}{(3+x)^5} and so we actually have the error between \log 2 and the approximation using Gauss quadrature with n=1 is

\begin{align} \log 2 - \left( \frac1{ 3 - \frac{\sqrt{3}}3 } + \frac1{ 3 + \frac{\sqrt{3}}3 } \right) \end{align}

and belongs to [\frac{24}{135} \min_{x\in [-1,1]} (3+x)^{-5}, \frac{24}{135} \max_{x\in [-1,1]} (3+x)^{-5}] = [ \frac1{5760}, \frac{1}{180} ] \supset [ 1.74 \times 10^{-4}, 5.55 \times 10^{-3} ]. The actual error belongs to this interval!

Exercise.

Compare this to rectangular and its composite rule with 2 sub-intervals and the trapezoid and Simpson rules.

15.5 Error estimates

If w_j \geq 0 and the quadrature rule is exact for all polynomials of degree N (this is true for Composite Trapezoid rule with N = 1 and Gaussian quadrature with N = 2n+1), then

\begin{align} &\left| \int_{-1}^1 f - \sum_{j=0}^n w_j f(x_j) \right| \\ % &= \min_{q \in \mathcal P_{N}}\left| \int_{-1}^1 (f - q) - \sum_{j=0}^n w_j \left( f(x_j) - q(x_j)\right) \right| \\ % &\leq 4 \min_{q \in \mathcal P_{N}} \left\| f - q \right\|_{L^\infty([-1,1])} \end{align}

c.f. Approximation theory [we will come back this this].

15.6 How to compute the roots of P_{n+1}?

Recall, P_0(x) = 1, P_1(x) = x, and

\begin{align} xP_n(x) = a_n P_n(x) + P_{n+1}(x) + b_n P_{n-1}(x) \end{align}

It turns out that that X = \{ \text{zeros of } P_{n+1} \} is the set of eigenvalues of a tridiagonal matrix T (see: Assignment 7). Next chapter, we will see methods for finding the eigenvalues of this matrix!

15.7 Exercises

Exercise 1.

Show that Gauss quadrature with nodes X = \{x_0, \dots, x_n\} is never exact for all polynomials of degree N > 2n+1.

Hint: Apply the quadrature rule to \ell_X(x)^2.

Exercise 2.

  1. Write down the Gaussian quadrature rule with n=2.
  2. What is the absolute error in approximating \log 2 using the n=2 scheme and following the same arguments as in class?
  3. How would you prove an error bound for the n=2 rule in general?

Exercise 3.

Try to understand the proof that roots of Legendre polynomials are contained in [-1,1].

15.8 Next…

  • Assignment 7: I’ll set A7 on Gaussian quadrature for you to do in the week beginning November 10
  • Next week we start Chapter 5: Solving Linear Equations - it would be very useful if you read through Section 6.1 of Burden, Faires, and Burden (2015) before the start of class on Monday.
  • Presentations: I will make a list of ideas for topics that you can present on. You’ll then be able to sign up to give a presentation. If you have a topic in mind, let me know and I can advise if your topic is suitable.
Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.